PROGRAM Main
  USE commonvars
  USE simulated_anneal
  USE Nelder_Mead
  USE Random

  IMPLICIT NONE
  INCLUDE 'mpif.h'

  ! SIMULATED ANNEALING
  REAL(8)            :: param(nparam), pp(nparam), fval, t1, VarCov(nparam,nparam)
  INTEGER, PARAMETER :: neps = 3
  REAL(8)            :: xopt(nparam), c(nparam) = 2d0, vm(nparam), t = 5.0, eps = 0.05, rt = 0.85, fopt
  INTEGER            :: ns = 20, nt = 30, nfcnev, iseed1 = 50+2, iseed2 = 80+2, maxevl = 100000, nacc, nobds, iiprint = 3, ier
  LOGICAL            :: ismax = .FALSE.
  REAL(8)            :: s, ttt, high_wtp, ch_edu, mu_finet, mag_ch_vec(nmag)
  INTEGER            :: i, j, tt, k, ii, jj, kk, Np_sim = 500, Nvar_sim, iboot, ntasks, nprocs_per_task, Nboot_tot, Npolicies, Nwageinc, Nwage_tries, nprau, n_chfine, n_prconv
  REAL(8), ALLOCATABLE :: simoutt(:,:,:,:), param_draws(:,:), pr_audit_vec(:), pc_cost_prau_vec(:), pc_cost_prau_3rd_vec(:), n_times_wage_vec(:), pc_cost_wage_vec(:,:), &
                          per_ch_fine_vec(:), prob_conv_vec(:)
  ! Nelder-Mead
  REAL(8)            :: object, simp, step(nparam), stopcr, var(nparam), temp, temp_vec(6), temp_age_vec(2)
  INTEGER            :: iiiprint, iquad, maxf, nloop, i_opt(6)
  LOGICAL            :: first
  EXTERNAL           :: objfcn, read_ini_data, fcn, standarderrors, standarderrors_subset, standarderrors_loop, set_policy_vars

  !Initialize the MPI execution environment
  CALL MPI_INIT(ierr)
  !Determines the size of the group associated with a communictor
  CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nprocs, ierr)
  !Determines the rank of each processor
  CALL MPI_COMM_RANK(MPI_COMM_WORLD, myrank, ierr)

  IF (myrank == 0) write(*,*) 'no. of processors used:', nprocs

  ALLOCATE(simoutt(Nactmun,Nsim,Nsimterms,Nvar))
  simoutt = 0d0

  ALLOCATE(av_qcons_vec(Tend,Nmun,Npaudits))
  av_qcons_vec = 0d0
  
  DEALLOCATE(simoutt)
  
  DO i = 1, Nwealth
     IF (wealth(i) == 0.0d0) THEN
        w0 = i
        EXIT
     END IF
  END DO

!!$  stealing = fr_stoleni*funds(CEILING(DBLE(Nfunds)/2d0))
  stealing = fr_stoleni

!!$  qconsvec(1,:) = qconsvec1-LOG(popsize(1))
!!$  qconsvec(2,:) = qconsvec2-LOG(popsize(2))
!!$  qconsvec(3,:) = qconsvec3-LOG(popsize(3))
  qconsvec(1,:) = qconsvec1
  qconsvec(2,:) = qconsvec2
  qconsvec(3,:) = qconsvec3
     
  Nnterms = 2
  n_times_wage = 1d0
  flag_no_run = 0
  paudit(0:1,1:Npaudits) = RESHAPE((/ 0.95d0, 0.05d0, 0.832d0, 0.168d0 /), SHAPE = (/2, Npaudits/)) ! (/ prob no audit, prob audit /)
  flag_audit_nterm = 0
 
  wage_mayor(:) = n_times_wage*(/ 3233.379d0, 4267.472d0, 5076.549d0 /)/norm
  paudit_nterm(0:1) = (/ 0.832d0, 0.168d0 /)
  
  ! Compute the machine epsilon used in probit
  s = 1d0
  do k=1,100
     s = 0.5d0*s
     ttt = s + 1.0d0

     if (ttt .le. 1.0d0) then
        s = 2.0d0*s
        ttt = 1.0d0/2.0d0**(k-1)
        IF (myrank == 0) write(*,*) 'machine precision k-1,s,ttt', k-1,s,ttt
        machine_eps = s
        exit
     endif
  end do

  IF (myrank == 0) write(*,*) 'machine_eps', machine_eps

  DO i = 1, Npaudits
     IF (myrank == 0) write(*,*) 'p. no aud, p. aud, i', paudit(0:1,i)
  END DO
     IF (myrank == 0) write(*,*) 'paudit_nterm', paudit_nterm(0:1)

  IF (myrank == 0) write(*,*) 'flag_sim', flag_sim
  
  ALLOCATE(udraw(Nactmun, Nsim, 1:Nsimterms, 6), age_udraw(Nactmun, Nsim, 1:Nsimterms, 2))
  ALLOCATE(fine_norm(Nactmun, Nsim, 1:Nsimterms), ability_norm(Nactmun, Nsim, 1:Nsimterms), theta_norm(Nactmun, Nsim, 1:Nsimterms))
  ALLOCATE(rnorm(Nactmun, Nsim, 2:Nsimterms))
  ALLOCATE(rnormprob(Nactmun, Nsim, 1:Nsimterms))
  ALLOCATE(zzpr_shock(Nactmun, Nsim, 1:Nsimterms))
  ALLOCATE(subutil(Tend))

  udraw = 0d0
  age_udraw = 0d0
  fine_norm = 0d0
  theta_norm = 0d0
  rnorm = 0d0
  rnormprob = 0d0
    
  !Generate draws from a random normal for the simulated data
  IF (flag_drawrandom == 1) THEN

     DO i = 1, Nactmun
        DO j = 1, Nsim
           DO tt = 2, Nsimterms
              
              !In Sample_Normal the second argument must be the variance
              rnorm(i,j,tt) = Sample_Normal(0.0d0,1.0d0)
              
           END DO
        END DO
     END DO
     
     DO i = 1, Nactmun
        DO j = 1, Nsim
           DO tt = 1, Nsimterms
              
              rnormprob(i,j,tt) = Sample_Normal(0.0d0,1.0d0)
              
           END DO
        END DO
     END DO

     DO i = 1, Nactmun
        DO j = 1, Nsim
           DO tt = 1, Nsimterms

              !In Sample_Normal the second argument must be the variance
              zzpr_shock(i,j,tt) = Sample_Normal(0.0d0,1.0d0)

           END DO
        END DO
     END DO

     !Generate the random uniform shocks for the simulated data
     DO i = 1, Nactmun
        DO j = 1, Nsim
           DO tt = 1, Nsimterms
              DO ii = 1, 6
                 udraw(i,j,tt,ii) = Sample_Uniform(0d0,1d0)
              END DO
           END DO
        END DO
     END DO

     !Generate the random uniform shocks for the simulated data
     DO i = 1, Nactmun
        DO j = 1, Nsim
           DO tt = 1, Nsimterms
              DO ii = 1, 2
                 age_udraw(i,j,tt,ii) = Sample_Uniform(0d0,1d0)
              END DO
           END DO
        END DO
     END DO
     
     DO i = 1, Nactmun
        DO j = 1, Nsim
           DO tt = 1, Nsimterms

              !In Sample_Normal the second argument must be the variance
              fine_norm(i,j,tt) = Sample_Normal(0d0,1d0)

           END DO
        END DO
     END DO

     DO i = 1, Nactmun
        DO j = 1, Nsim
           DO tt = 1, Nsimterms

              !In Sample_Normal the second argument must be the variance
              ability_norm(i,j,tt) = Sample_Normal(0d0,1d0)

           END DO
        END DO
     END DO

     DO i = 1, Nactmun
        DO j = 1, Nsim
           DO tt = 1, Nsimterms
              
              !In Sample_Normal the second argument must be the variance
              theta_norm(i,j,tt) = Sample_Normal(0d0,1d0)

           END DO
        END DO
     END DO

!!$     !Generate the random uniform shocks for the simulated data
!!$     DO i = 1, Nactmun
!!$        DO j = 1, Nboot_max
!!$           uboot(j,i) = Sample_Uniform(0d0,1d0)
!!$        END DO
!!$     END DO
     
     IF (myrank == 0) THEN

        IF (Nsim < 101) THEN

           OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/mayor_rnorm.dat")
           DO i = 1, Nactmun
              DO j = 1, Nsim
                 DO tt = 2, Nsimterms
                    
                    WRITE(41,*) rnorm(i,j,tt)
                    
                 ENDDO
              END DO
           END DO
           CLOSE(41)
           
           OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/mayor_rnormprob.dat")
           DO i = 1, Nactmun
              DO j = 1, Nsim
                 DO tt = 1, Nsimterms
                    
                    WRITE(41,*) rnormprob(i,j,tt)
                    
                 ENDDO
              END DO
           END DO
           CLOSE(41)
           
           OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/zzpr_shocks.dat")
           DO i = 1, Nactmun
              DO j = 1, Nsim
                 DO tt = 1, Nsimterms
                    
                    WRITE(41,*) zzpr_shock(i,j,tt)
                    
                 ENDDO
              END DO
           END DO
           CLOSE(41)
           
           OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/mayor_udraw.dat")
           DO i = 1, Nactmun
              DO j = 1, Nsim
                 DO tt = 1, Nsimterms
                    
                    WRITE(41,*) udraw(i,j,tt,:)
                    
                 ENDDO
              END DO
           END DO
           CLOSE(41)

           OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/age_udraw.dat")
           DO i = 1, Nactmun
              DO j = 1, Nsim
                 DO tt = 1, Nsimterms
                    
                    WRITE(41,*) age_udraw(i,j,tt,:)
                    
                 ENDDO
              END DO
           END DO
           CLOSE(41)

           OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/fine_norm.dat")
           DO i = 1, Nactmun
              DO j = 1, Nsim
                 DO tt = 1, Nsimterms
                    WRITE(41,*) fine_norm(i,j,tt)
                 END DO
              END DO
           END DO
           CLOSE(41)
           
           OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/ability_norm.dat")
           DO i = 1, Nactmun
              DO j = 1, Nsim
                 DO tt = 1, Nsimterms
                    WRITE(41,*) ability_norm(i,j,tt)
                 END DO
              END DO
           END DO
           CLOSE(41)

           OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/theta_norm.dat")
           DO i = 1, Nactmun
              DO j = 1, Nsim
                 DO tt = 1, Nsimterms
                    WRITE(41,*) theta_norm(i,j,tt)
                 END DO
              END DO
           END DO
           CLOSE(41)

        ELSE
           
           OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/mayor_rnorm_500.dat")
           DO i = 1, Nactmun
              DO j = 1, Nsim
                 DO tt = 2, Nsimterms
                    
                    WRITE(41,*) rnorm(i,j,tt)
                    
                 ENDDO
              END DO
           END DO
           CLOSE(41)
           
           OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/mayor_rnormprob_500.dat")
           DO i = 1, Nactmun
              DO j = 1, Nsim
                 DO tt = 1, Nsimterms
                    
                    WRITE(41,*) rnormprob(i,j,tt)
                    
                 ENDDO
              END DO
           END DO
           CLOSE(41)
           
           OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/zzpr_shocks_500.dat")
           DO i = 1, Nactmun
              DO j = 1, Nsim
                 DO tt = 1, Nsimterms
                    
                    WRITE(41,*) zzpr_shock(i,j,tt)
                    
                 ENDDO
              END DO
           END DO
           CLOSE(41)
           
           OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/mayor_udraw_500.dat")
           DO i = 1, Nactmun
              DO j = 1, Nsim
                 DO tt = 1, Nsimterms
                    
                    WRITE(41,*) udraw(i,j,tt,:)
                    
                 ENDDO
              END DO
           END DO
           CLOSE(41)
           
           OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/fine_norm_500.dat")
           DO i = 1, Nactmun
              DO j = 1, Nsim
                 DO tt = 1, Nsimterms
                    WRITE(41,*) fine_norm(i,j,tt)
                 END DO
              END DO
           END DO
           CLOSE(41)
           
           OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/ability_norm_500.dat")
           DO i = 1, Nactmun
              DO j = 1, Nsim
                 DO tt = 1, Nsimterms
                    WRITE(41,*) ability_norm(i,j,tt)
                 END DO
              END DO
           END DO
           CLOSE(41)

           OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/theta_norm_500.dat")
           DO i = 1, Nactmun
              DO j = 1, Nsim
                 DO tt = 1, Nsimterms
                    WRITE(41,*) theta_norm(i,j,tt)
                 END DO
              END DO
           END DO
           CLOSE(41)

        END IF

!!$        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/random_bootstrap_numbers.dat")
!!$        DO i = 1, Nactmun
!!$           DO j = 1, Nboot_max
!!$              WRITE(41,*) uboot(j,i)
!!$           END DO
!!$        END DO
!!$        CLOSE(41)
        
     END IF

  ELSE
     
     IF (Nsim < 101) THEN
        
        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/mayor_rnorm.dat")
        DO i = 1, Nactmun
           DO j = 1, Nsim
              DO tt = 2, Nsimterms
                 
                 READ(41,*) temp
                 rnorm(i,j,tt) = temp
                 
              ENDDO
           END DO
        END DO
        CLOSE(41)
        
        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/mayor_rnormprob.dat")
        DO i = 1, Nactmun
           DO j = 1, Nsim
              DO tt = 1, Nsimterms
                 
                 READ(41,*) temp
                 rnormprob(i,j,tt) = temp
                 
              ENDDO
           END DO
        END DO
        CLOSE(41)
        
        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/zzpr_shocks.dat")
        DO i = 1, Nactmun
           DO j = 1, Nsim
              DO tt = 1, Nsimterms
                 
                 READ(41,*) temp
                 zzpr_shock(i,j,tt) = temp
                 
              ENDDO
           END DO
        END DO
        CLOSE(41)
        
        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/mayor_udraw.dat")
        DO i = 1, Nactmun
           DO j = 1, Nsim
              DO tt = 1, Nsimterms
                 
                 READ(41,*) temp_vec
                 udraw(i,j,tt,:) = temp_vec(:)
                 
              ENDDO
           END DO
        END DO
        CLOSE(41)

        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/age_udraw.dat")
        DO i = 1, Nactmun
           DO j = 1, Nsim
              DO tt = 1, Nsimterms
                 
                 READ(41,*) temp_age_vec
                 age_udraw(i,j,tt,:) = temp_age_vec(:)
                 
              ENDDO
           END DO
        END DO
        CLOSE(41)

        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/fine_norm.dat")
        DO i = 1, Nactmun
           DO j = 1, Nsim
              DO tt = 1, Nsimterms
                 READ(41,*) temp
                 fine_norm(i,j,tt) = temp
              END DO
           END DO
        END DO
        CLOSE(41)
        
        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/ability_norm.dat")
        DO i = 1, Nactmun
           DO j = 1, Nsim
              DO tt = 1, Nsimterms
                 READ(41,*) temp
                 ability_norm(i,j,tt) = temp
              END DO
           END DO
        END DO
        CLOSE(41)

        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/theta_norm.dat")
        DO i = 1, Nactmun
           DO j = 1, Nsim
              DO tt = 1, Nsimterms
                 READ(41,*) temp
                 theta_norm(i,j,tt) = temp
              END DO
           END DO
        END DO
        CLOSE(41)

     ELSE

        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/mayor_rnorm_500.dat")
        DO i = 1, Nactmun
           DO j = 1, Nsim
              DO tt = 2, Nsimterms
                 
                 READ(41,*) temp
                 rnorm(i,j,tt) = temp
                 
              ENDDO
           END DO
        END DO
        CLOSE(41)
        
        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/mayor_rnormprob_500.dat")
        DO i = 1, Nactmun
           DO j = 1, Nsim
              DO tt = 1, Nsimterms
                 
                 READ(41,*) temp
                 rnormprob(i,j,tt) = temp
                 
              ENDDO
           END DO
        END DO
        CLOSE(41)
        
        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/zzpr_shocks_500.dat")
        DO i = 1, Nactmun
           DO j = 1, Nsim
              DO tt = 1, Nsimterms
                 
                 READ(41,*) temp
                 zzpr_shock(i,j,tt) = temp
                 
              ENDDO
           END DO
        END DO
        CLOSE(41)
        
        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/mayor_udraw_500.dat")
        DO i = 1, Nactmun
           DO j = 1, Nsim
              DO tt = 1, Nsimterms
                 
                 READ(41,*) temp_vec
                 udraw(i,j,tt,:) = temp_vec(:)
                 
              ENDDO
           END DO
        END DO
        CLOSE(41)
        
        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/fine_norm_500.dat")
        DO i = 1, Nactmun
           DO j = 1, Nsim
              DO tt = 1, Nsimterms
                 READ(41,*) temp
                 fine_norm(i,j,tt) = temp
              END DO
           END DO
        END DO
        CLOSE(41)
        
        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/ability_norm_500.dat")
        DO i = 1, Nactmun
           DO j = 1, Nsim
              DO tt = 1, Nsimterms
                 READ(41,*) temp
                 ability_norm(i,j,tt) = temp
              END DO
           END DO
        END DO
        CLOSE(41)

        OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/theta_norm_500.dat")
        DO i = 1, Nactmun
           DO j = 1, Nsim
              DO tt = 1, Nsimterms
                 READ(41,*) temp
                 theta_norm(i,j,tt) = temp
              END DO
           END DO
        END DO
        CLOSE(41)
        
     END IF

!!$     OPEN(unit=41, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/random_bootstrap_numbers.dat")
!!$     DO i = 1, Nactmun
!!$        DO j = 1, Nboot_max
!!$           READ(41,*) temp
!!$           uboot(j,i) = temp
!!$        END DO
!!$     END DO
!!$     CLOSE(41)

  END IF
  
  min_ability = 0d0 !MINVAL(ability_norm(:,:,:))
  
  !Compute the size of all loop
  IF (flag_policy_evaluation == 0) THEN
  
     ALLOCATE(ubp(nparam), lbp(nparam))
     
     OPEN(22, file='inputs.txt')
     READ(22,*) pp
     READ(22,*) ubp
     READ(22,*) lbp
     READ(22,*) factor_param
     CLOSE(22)
     
     DO i = 1, nparam
        IF (myrank == 0) write(*,'(a,1i4,2f14.7)') 'i, pp, factor_param', i, pp(i), factor_param(i)
     END DO
     
     DO i = 1, nparam
        IF (myrank == 0) write(*,'(a,1i4,3f20.10)') 'i, lbp, pp, ubp', i, lbp(i), pp(i), ubp(i)
     END DO
     
     DO i = 1, nparam
        param(i) = LOG((pp(i)-lbp(i))/(ubp(i)-pp(i)))
        IF (myrank == 0) write(*,'(a,1i4,1f14.7)') 'i, param', i, param(i)
     END DO
     
     DO i = 1, nparam
        vm(i) = 10d0
     END DO

     CALL set_policy_vars
          
     !Read in the matrix of initial conditions
     !the matrix characterizing the nonparameteric distribution
     !the nonparametric distribution for emax
     CALL read_ini_data

     IF (flag_simulate_data == 0 .AND. flag_compstderrs == 0) THEN
        
        IF (flag_sa == 1) THEN
           
           !-----------------------------------------------------------------
           ! estimation by SA
           !-----------------------------------------------------------------
           
           param = pp
           
           DO i = 1, nparam
              if (myrank == 0) write(*,*) 'u,l,param', ubp(i), lbp(i), param(i)
           END DO
           
           IF (myrank == 0) THEN
              write(*,*) 'param', param
           END IF
           
           IF (flag_sim == 1) THEN
              
              flag_sim_data = 1
              CALL fcn(nparam, param, fval)
              flag_sim_data = 0
              
              STOP
           END IF
           
           ! Estimation by simulated annealing
           CALL sa(nparam, param, ismax, rt, eps, ns, nt, neps, maxevl, lbp, ubp, c, iiprint, iseed1,  &
                iseed2, t, vm, xopt, fopt, nacc, nfcnev, nobds, ier)
           
           DEALLOCATE(ubp,lbp)
           
        ELSE
           
           !------------------------------------------------------------------
           ! estimation by Nelder-Mead
           !-----------------------------------------------------------------
           
           ! step sizes
           DO i = 1, nparam
              step(i) = 1d0
           END DO
           
           ! max # of function evaluations
           maxf = 10000
           ! print every iteraton
           iiiprint = 1
           ! stopping criterion
           stopcr = 1d-5 !1d-4
           nloop = 2*nparam
           ! don't fit a quadratic surface
           iquad = 0
           ! accuracy 9 digits
           simp = 1d-6
           ! declare this is first time
           first = .true.
           
           CALL minim(param, step, nparam, object, maxf, iiiprint, stopcr, nloop, iquad, simp, var, objfcn, ier)
           
        END IF
        
     ELSE IF (flag_compstderrs == 1) THEN

        CALL standarderrors(nparam,param(:))

     ELSE IF (flag_simulate_data == 1) THEN

        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !We then run the code to generate the data we are interested in!
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

        flag_estimation = 1
        CALL set_policy_vars

        ! We run the code
        CALL fcn(nparam, param, fval)

        flag_estimation = 0
        DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

        flag_base_policy = 1
        CALL set_policy_vars

        ! We run the code
        CALL fcn(nparam, param, fval)

        flag_base_policy = 0
        DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

        flag_1term_policy = 1
        CALL set_policy_vars

        ! We run the code
        CALL fcn(nparam, param, fval)

        flag_1term_policy = 0
        DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

        flag_base_policy = 1; flag_no_selec_sav = 1
        CALL set_policy_vars

        ! We run the code
        CALL fcn(nparam, param, fval)

        flag_base_policy = 0; flag_no_selec_sav = 0
        DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

        flag_1term_policy = 1; flag_no_selec_mun = 1; flag_no_selec_funds = 1; flag_no_selec_zzpr = 1; flag_no_selec_ab = 1; flag_no_selec_educ = 1; flag_no_selec_age = 1
        CALL set_policy_vars

        ! We run the code
        CALL fcn(nparam, param, fval)

        flag_1term_policy = 0; flag_no_selec_mun = 0; flag_no_selec_funds = 0; flag_no_selec_zzpr = 0; flag_no_selec_ab = 0; flag_no_selec_educ = 0; flag_no_selec_age = 0
        DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

        !We then run the code for the policy we are interested in
        flag_norun_policy = 1
        CALL set_policy_vars
        
        ! We run the policies
        CALL fcn(nparam, param, fval)
        
        flag_norun_policy = 0
        DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

        !We then run the code for the policy we are interested in
        flag_terms_policy = 1
        CALL set_policy_vars
        
        ! We run the policies
        CALL fcn(nparam, param, fval)

        flag_terms_policy = 0
        DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

     END IF

  ELSE IF (flag_policy_evaluation == 1 .AND. flag_optimal_auditprob == 1) THEN

     ! cost_policy = 1 when we account for the cost of audits or increased wages
     cost_policy = 0
     ! Nvar_sim is equal to the number of variable we want to save independently of terms (5) plus the highest number of terms we have in all our policies (3) 
     Nvar_sim = 5+3
     ! Number of audit probabilities we consider
     Nprau = 18
     ! Number of wage increases we consider
     Nwageinc = 9
     n_chfine = 20
     n_prconv = 3
     ! Npolicies = 10 when we simulate the policies
     ! Npolicies = 3*Nprau+Nwage_tries
     Npolicies = 7*Nprau
     ! Npolicies = Nwageinc

     ! For this part of the program we don't neet to create groups of processors
     Nboot_tot = 1
     iboot = 1
     nboot = 1
     first_proc = 1
     
     ALLOCATE(bootstrap_vect(Npolicies,Nboot_tot,Nvar_sim), bootstrap_vec(Npolicies,Nboot_tot,Nvar_sim))
     bootstrap_vect = 0d0
     bootstrap_vec = 0d0
     ALLOCATE(notmayor_value_t(Nactmun,Nsim,Nsimterms), notmayor_value(Nactmun,Nsim,Nsimterms), ELNqcons_base_vec_t(Nactmun,Nsim,Nsimterms), ELNqcons_base_vec(Nactmun,Nsim,Nsimterms))
     notmayor_value_t = 0d0
     notmayor_value = 0d0

     ALLOCATE(pr_audit_vec(nprau), pc_cost_prau_vec(nprau), pc_cost_prau_3rd_vec(nprau), n_times_wage_vec(nwageinc), pc_cost_wage_vec(nmun,nwageinc))
     pr_audit_vec(:) = (/ ( 0.02d0+0.01d0*DBLE(i) , i=1,nprau) /)
     ! Total Cost of the Audit Policy
     pc_cost_prau_vec(:) = (/ 0.38d0, 0.5d0, 0.63d0, 0.75d0, 0.88d0, 1.01d0, 1.13d0, 1.26d0, 1.38d0, 1.51d0, 1.63d0, 1.76d0, 1.89d0, 2.01d0, 2.14d0, 2.26d0, 2.39d0, 2.52d0 /)
     ! Increase in Cost of the Audit policy (0.63 is the cost if the pr_audit is 0.5)
     pc_cost_prau_vec = pc_cost_prau_vec - 0.63d0
     ! Cost per tax-payer of an audit policy that changes the probability of audits only in the third term relative to the base case (5% audit all terms)
     ! We subtract 0.02898d0, cost of the 3rd term audit policy if 3rd term mayors are audit with a 5% probability; we therefore have the increase relative to the base case
     ! In the vector below, We start with the cost for the audit prob 82% because the optimal prob is 85% and we need a vector with 18 elements to be consistent with above
     pc_cost_prau_3rd_vec(:) = (/ 0.475640d0, 0.480240d0, 0.487130d0, 0.491730d0, 0.498620d0, 0.503220d0, 0.510110d0, 0.514710d0, 0.521600d0, 0.526190d0, .533090d0, 0.537680d0, 0.544580d0, 0.549170d0, 0.556070d0, 0.560660d0, 0.567560d0, 0.572150d0 /) - 0.02898d0

     ! To compute the optimal wage policy
!!$     n_times_wage_vec(:) = (/ 1d0, 1.5d0, 2d0, 2.5d0, 3d0, 3.5d0, 4d0, 4.5d0, 5d0 /)
     n_times_wage_vec(:) = (/ 1d0, 4.2d0, 4.3d0, 4.4d0, 4.5d0, 4.6d0, 4.7d0, 4.8d0, 5d0 /)
     pc_cost_wage_vec(1,:) =  0.7059904d0*(n_times_wage_vec(:) - 1d0)
     pc_cost_wage_vec(2,:) =  0.2450834d0*(n_times_wage_vec(:) - 1d0)
     pc_cost_wage_vec(3,:) =  0.0580153d0*(n_times_wage_vec(:) - 1d0)
     
     ! WE COMPUTE THE OPTIMAL AUDIT PROBABILITY FOR DIFFERENT POLICIES

     !Read in the matrix of initial conditions
     !the matrix characterizing the nonparameteric distribution
     !the nonparametric distribution for emax
     CALL read_ini_data

     ALLOCATE(ubp(nparam), lbp(nparam))

     OPEN(22, file='inputs.txt')
     READ(22,*) pp
     READ(22,*) ubp
     READ(22,*) lbp
     READ(22,*) factor_param
     CLOSE(22)
     
     DO i = 1, nparam
        IF (myrank == 0) write(*,'(a,1i4,2f14.7)') 'i, pp, factor_param', i, pp(i), factor_param(i)
     END DO
     
     DO i = 1, nparam
        IF (myrank == 1) write(*,'(a,1i4,3f20.10)') 'i, lbp, pp, ubp', i, lbp(i), pp(i), ubp(i)
     END DO
     
     DO i = 1, nparam
        param(i) = LOG((pp(i)-lbp(i))/(ubp(i)-pp(i)))
     END DO
     
     ! We need param = pp and flag_sa = 1 to use fcn, if not the parameters are transformed
     param = pp
     flag_sa = 1

     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !We first run the code for the BASE CASE to get the base past mayor values!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
     flag_base_policy = 1
     ELNqcons_base_vec_t = 0d0
     ELNqcons_base_vec = 0d0
     policy_pos = Npolicies
     CALL set_policy_vars
     
     ! We run the policies
     CALL fcn(nparam, param, fval)
     
     flag_base_policy = 0
     DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !We then run the code for the OPTIMAL AUDIT POLICY!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     ! We compute the optimal audit policy
     flag_optimal_prau_policy = 1
     cost_policy = 1
     DO i = 1, nprau

        pr_audit_out = pr_audit_vec(i)
        pc_cost_prau_out = pc_cost_prau_vec(i)
        !We then run the code for the policy we are interested in
        policy_pos = i
        CALL set_policy_vars
        
        ! We run the policies
        CALL fcn(nparam, param, fval)
        
        DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)
        
     END DO

     IF (myrank == 0) write(*,'(a)') 'END OPTIMAL AUDIT POLICY'

     flag_optimal_prau_policy = 0
     cost_policy = 0
     
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !We then run the code for the TERM LIMITS WITH OPTIMAL AUDIT POLICY!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     jj = nprau
     ! We compute the optimal audit policy
     flag_term_optimal_prau_policy = 1
     cost_policy = 1
     DO i = 1, nprau
        
        pr_audit_out = pr_audit_vec(i)
        pc_cost_prau_out = pc_cost_prau_vec(i)
        !We then run the code for the policy we are interested in
        policy_pos = jj + i
        CALL set_policy_vars
        
        ! We run the policies
        CALL fcn(nparam, param, fval)
        
        DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)
        
     END DO

     IF (myrank == 0) write(*,'(a)') 'END TERM LIMITS WITH OPTIMAL AUDIT POLICY'

     flag_term_optimal_prau_policy = 0
     cost_policy = 0
     
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !We then run the code for the NO-RUN WITH OPTIMAL AUDIT POLICY!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     jj = 2*nprau 
     ! We compute the optimal audit policy
     flag_norun_optimal_prau_policy = 1
     cost_policy = 1
     DO i = 1, nprau
        
        pr_audit_out = pr_audit_vec(i)
        pc_cost_prau_out = pc_cost_prau_vec(i)
        !We then run the code for the policy we are interested in
        policy_pos = jj + i
        CALL set_policy_vars
        
        ! We run the policies
        CALL fcn(nparam, param, fval)
        
        DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)
        
     END DO

     IF (myrank == 0) write(*,'(a)') 'END NO-RUN WITH OPTIMAL AUDIT POLICY'

     flag_norun_optimal_prau_policy = 0
     cost_policy = 0
     
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !We then run the code for the TERM AND NO-RUN WITH OPTIMAL AUDIT POLICY!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     jj = 3*nprau 
     ! We compute the optimal audit policy
     flag_term_norun_optimal_prau_policy = 1
     cost_policy = 1
     DO i = 1, nprau
        
        pr_audit_out = pr_audit_vec(i)
        pc_cost_prau_out = pc_cost_prau_vec(i)
        !We then run the code for the policy we are interested in
        policy_pos = jj + i
        CALL set_policy_vars
        
        ! We run the policies
        CALL fcn(nparam, param, fval)
        
        DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)
        
     END DO
     
     IF (myrank == 0) write(*,'(a)') 'END TERM AND NO-RUN WITH OPTIMAL AUDIT POLICY'

     flag_term_norun_optimal_prau_policy = 0
     cost_policy = 0           

     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !We then run the code for the TERM AN NO-RUN WITH OPTIMAL AUDIT POLICY ONLY IN THE LAST TERM!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     jj = 4*nprau 
     ! We compute the optimal audit policy
     flag_term_norun_optimal_prau_lastterm_policy = 1
     cost_policy = 1
     DO i = 1, nprau

        ! We have 0.8+pr_audit_vec(i) to have as the starting point an audit prob = 83% (the optimal is 85%)
        pr_audit_out = 0.79+pr_audit_vec(i)
        !pr_audit_out = pr_audit_vec(i)
        pc_cost_prau_out = pc_cost_prau_3rd_vec(i)
        !We then run the code for the policy we are interested in
        policy_pos = jj + i
        CALL set_policy_vars
        
        ! We run the policies
        CALL fcn(nparam, param, fval)
        
        DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)
        
     END DO
     
     IF (myrank == 0) write(*,'(a)') 'END TERM AND NO-RUN WITH OPTIMAL AUDIT POLICY ONLY IN THE LAST TERM!'

     flag_term_norun_optimal_prau_lastterm_policy = 0
     cost_policy = 0           

     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !We then run the code for the WAGE POLICY WITH OPTIMAL AUDIT PROBABILITY!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     jj = 5*nprau
     ! We compute the optimal audit policy
     flag_wage_optimal_prau_policy = 1
     cost_policy = 1
     
     DO i = 1, nprau
        
        pr_audit_out = pr_audit_vec(i)
        pc_cost_prau_out = pc_cost_prau_vec(i)
        
        !We then run the code for the policy we are interested in
        policy_pos = jj + i
        CALL set_policy_vars
        
        ! We run the policies
        CALL fcn(nparam, param, fval)
        
        DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)
        
     END DO

     IF (myrank == 0) write(*,'(a)') 'END WAGE POLICY WITH OPTIMAL AUDIT PROBABILITY'

     flag_wage_optimal_prau_policy = 0
     cost_policy = 0

     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !We then run the code for the OPTIMAL WAGE POLICY !
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     jj = 6*nprau
!!$     jj = 0
     ! We compute the optimal audit policy
     flag_optimal_wage_policy = 1
     cost_policy = 1
     
     DO i = 1, nwageinc
        
        n_times_wage_out = n_times_wage_vec(i)
        pc_cost_wage_out(:) = pc_cost_wage_vec(:,i)
        
        !We then run the code for the policy we are interested in
        policy_pos = jj + i
        CALL set_policy_vars
        
        ! We run the policies
        CALL fcn(nparam, param, fval)
        
        DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)
        
     END DO

     IF (myrank == 0) write(*,'(a)') 'END OPTIMAL WAGE POLICY'

     flag_optimal_wage_policy = 0
     cost_policy = 0

!!$     CALL MPI_ALLREDUCE(bootstrap_vect(1,1,1), bootstrap_vec(1,1,1), Npolicies*Nboot_tot*Nvar_sim, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
     bootstrap_vec = bootstrap_vect
     
     IF (myrank == 0) THEN
        
        high_wtp = -1000d0
        DO i = 1, nprau
           IF (bootstrap_vec(i,1,4) > high_wtp) THEN
              high_wtp = bootstrap_vec(i,1,4)
              i_opt(1) = i
              write(*,'(a,1i4,10f16.10)') 'HIGH WTP Audit policy', i, pr_audit_vec(i), pc_cost_prau_vec(i), bootstrap_vec(i,1,:)
              !write(*,'(a,1i4,9f20.10)') 'HIGH WTP policy', i, ch_mu_a, bootstrap_vec(i,1,:)
           ELSE
              write(*,'(a,1i4,10f16.10)') 'Audit policy', i, pr_audit_vec(i), pc_cost_prau_vec(i), bootstrap_vec(i,1,:)
              !write(*,'(a,1i4,9f20.10)') 'policy', i, ch_mu_a, bootstrap_vec(i,1,:)
           END IF
        END DO
        high_wtp = -1000d0
        DO i = 1, nprau
           IF (bootstrap_vec(nprau+i,1,4) > high_wtp) THEN
              high_wtp = bootstrap_vec(nprau+i,1,4)
              i_opt(2) = i
              write(*,'(a,1i4,10f16.10)') 'HIGH WTP Term Limit policy', nprau+i, pr_audit_vec(i), pc_cost_prau_vec(i), bootstrap_vec(nprau+i,1,:)
           ELSE
              write(*,'(a,1i4,10f16.10)') 'Term Limit policy', nprau+i, pr_audit_vec(i), pc_cost_prau_vec(i), bootstrap_vec(nprau+i,1,:)
           END IF
        END DO
        high_wtp = -1000d0
        DO i = 1, nprau
           IF (bootstrap_vec(2*nprau+i,1,4) > high_wtp) THEN
              high_wtp = bootstrap_vec(2*nprau+i,1,4)
              i_opt(3) = i
              write(*,'(a,1i4,10f16.10)') 'HIGH WTP No-run policy', 2*nprau+i, pr_audit_vec(i), pc_cost_prau_vec(i), bootstrap_vec(2*nprau+i,1,:)
           ELSE
              write(*,'(a,1i4,10f16.10)') ' No-run policy', 2*nprau+i, pr_audit_vec(i), pc_cost_prau_vec(i), bootstrap_vec(2*nprau+i,1,:)
           END IF
        END DO
        high_wtp = -1000d0
        DO i = 1, nprau
           IF (bootstrap_vec(3*nprau+i,1,4) > high_wtp) THEN
              high_wtp = bootstrap_vec(3*nprau+i,1,4)
              i_opt(4) = i
              write(*,'(a,1i4,10f16.10)') 'HIGH WTP Term and No-run policy', 3*nprau+i, pr_audit_vec(i), pc_cost_prau_vec(i), bootstrap_vec(3*nprau+i,1,:)
           ELSE
              write(*,'(a,1i4,10f16.10)') 'Term and No-run policy', 3*nprau+i, pr_audit_vec(i), pc_cost_prau_vec(i), bootstrap_vec(3*nprau+i,1,:)
           END IF
        END DO
        high_wtp = -1000d0
        DO i = 1, nprau
           IF (bootstrap_vec(4*nprau+i,1,4) > high_wtp) THEN
              high_wtp = bootstrap_vec(4*nprau+i,1,4)
              i_opt(5) = i
              write(*,'(a,1i4,10f16.10)') 'HIGH WTP Term and No-run Policy 3rd', 4*nprau+i, pr_audit_vec(i), pc_cost_prau_3rd_vec(i), bootstrap_vec(4*nprau+i,1,:)
           ELSE
              write(*,'(a,1i4,10f16.10)') 'Term and No-run Policy 3rd', 4*nprau+i, pr_audit_vec(i), pc_cost_prau_3rd_vec(i), bootstrap_vec(4*nprau+i,1,:)
           END IF
        END DO
        high_wtp = -1000d0
        DO i = 1, nprau
           IF (bootstrap_vec(5*nprau+i,1,4) > high_wtp) THEN
              high_wtp = bootstrap_vec(5*nprau+i,1,4)
              i_opt(6) = i
              write(*,'(a,1i4,10f16.10)') 'HIGH WTP Wage Optimal Audit policy', 5*nprau+i, pr_audit_vec(i), pc_cost_prau_vec(i), bootstrap_vec(5*nprau+i,1,:)
           ELSE
              write(*,'(a,1i4,10f16.10)') 'policy', 5*nprau+i, pr_audit_vec(i), pc_cost_prau_vec(i), bootstrap_vec(5*nprau+i,1,:)
           END IF
        END DO
        high_wtp = -1000d0
        DO i = 1, nprau
           IF (bootstrap_vec(6*nprau+i,1,4) > high_wtp) THEN
              high_wtp = bootstrap_vec(6*nprau+i,1,4)
              i_opt(6) = i
              write(*,'(a,1i4,10f16.10)') 'HIGH WTP Optimal Wage policy', 6*nprau+i, pr_audit_vec(i), pc_cost_prau_vec(i), bootstrap_vec(6*nprau+i,1,:)
           ELSE
              write(*,'(a,1i4,10f16.10)') 'policy', 6*nprau+i, pr_audit_vec(i), pc_cost_prau_vec(i), bootstrap_vec(6*nprau+i,1,:)
           END IF
        END DO

!!$        DO i = 1, nwageinc
!!$           IF (bootstrap_vec(i,1,4) > high_wtp) THEN
!!$              high_wtp = bootstrap_vec(i,1,4)
!!$              i_opt(6) = i
!!$              write(*,'(a,1i4,10f16.10)') 'HIGH WTP Wage policy', i, pr_audit_vec(i), pc_cost_prau_vec(i), bootstrap_vec(i,1,:)
!!$           ELSE
!!$              write(*,'(a,1i4,10f16.10)') 'policy', i, pr_audit_vec(i), pc_cost_prau_vec(i), bootstrap_vec(i,1,:)
!!$           END IF
!!$        END DO

     END IF

!!$     OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/optimal_audit_probs.dat", STATUS='REPLACE')
!!$     WRITE(28,*) i_opt
!!$     CLOSE(28)
     
  ELSE IF (flag_policy_evaluation == 1 .AND. flag_optimal_auditprob == 0) THEN
     
     ! Read in the estimated parameters and standard errors
     OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/estparam_stderrs_tot.dat")
     READ(28,'(11f25.18)') pp(:)
     IF (myrank == 0) write(*,'(a,11f20.10)') 'param', pp
     DO i = 1, nparam
        READ(28,'(11f25.18)') VarCov(i,:)
     END DO

     ! Draw or read in the parameter shocks from a normal(0,1)
     ALLOCATE(param_draws(nparam,Np_sim))
     IF (flag_drawparam == 1) THEN

        DO j = 1, Np_sim
           ! Draw from a joint normal with mean pp and varcov VarCov
           param_draws(:,j) = Sample_Multivariate_Normal(pp,VarCov)
!!$           param_draws(5,j) =  param_draws(5,j) * 0.1d0 
           IF (myrank == 0 .AND. j <= 11) write(*,'(a,11f20.10)') 'param_draws(:,j)', param_draws(:,j)
        END DO

        IF (myrank == 0) THEN
           
           OPEN(unit=28, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/param_draws.dat")

           DO j = 1, Np_sim
              WRITE(28,*) param_draws(:,j)                    
           END DO
           CLOSE(28)

        END IF
        
     ELSE

        OPEN(unit=28, FILE = "/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/param_draws.dat")

        DO j = 1, Np_sim
           READ(28,*) param_draws(:,j)                    
        END DO
        CLOSE(28)

     END IF

     Nboot_tot = Nboot_max-Nboot_min+1
     ! cost_policy = 1 when we account for the cost of audits or increased wages
     cost_policy = 0
     ! Nvar_sim is equal to the number of variable we want to save independently of terms (5) plus the highest number of terms we have in all our policies (3) 
     Nvar_sim = 5+3
     ! Number of audit probabilities we consider
     Nprau = 18
     ! Number of wage increases we consider
     Nwage_tries = 13 !8
     n_chfine = 20
     n_prconv = 3
     ! Npolicies = 10 when we simulate the policies
     ! Npolicies = 3*Nprau+Nwage_tries
     Npolicies = 12 !5*Nprau !n_chfine+n_prconv+1 !5*Nprau !Nprau*Nwage_tries !Nprau+Nprau*Nwage_tries ! 11 
     
     ALLOCATE(bootstrap_vect(Npolicies,Nboot_tot,Nvar_sim), bootstrap_vec(Npolicies,Nboot_tot,Nvar_sim))
     bootstrap_vect = 0d0
     bootstrap_vec = 0d0
     ALLOCATE(notmayor_value_t(Nactmun,Nsim,Nsimterms), notmayor_value(Nactmun,Nsim,Nsimterms), ELNqcons_base_vec_t(Nactmun,Nsim,Nsimterms), ELNqcons_base_vec(Nactmun,Nsim,Nsimterms))
     notmayor_value_t = 0d0
     notmayor_value = 0d0

     ALLOCATE(pr_audit_vec(nprau), pc_cost_prau_vec(nprau), pc_cost_prau_3rd_vec(nprau))
     pr_audit_vec(:) = (/ ( 0.02d0+0.01d0*DBLE(i) , i=1,nprau) /)
     ! Total Cost of the Audit Policy
!!$     pc_cost_prau_vec(:) = (/ 0.75d0, 0.88d0, 1.01d0, 1.13d0, 1.26d0, 1.38d0, 1.51d0, 1.63d0, 1.76d0, 1.89d0, 2.01d0, 2.14d0, 2.26d0, 2.39d0, 2.52d0 /)
!!$     pc_cost_prau_vec(:) = (/ 0.63d0, 0.75d0, 0.88d0, 1.01d0, 1.13d0, 1.26d0, 1.38d0, 1.51d0, 1.63d0, 1.76d0, 1.89d0 /)
!!$     pc_cost_prau_vec(:) = (/ 0d0, 0.13d0, 0.25d0, 0.38d0, 0.5d0, 0.63d0, 0.75d0, 0.88d0, 1.01d0, 1.13d0, 1.26d0 /)
     pc_cost_prau_vec(:) = (/ 0.38d0, 0.5d0, 0.63d0, 0.75d0, 0.88d0, 1.01d0, 1.13d0, 1.26d0, 1.38d0, 1.51d0, 1.63d0, 1.76d0, 1.89d0, 2.01d0, 2.14d0, 2.26d0, 2.39d0, 2.52d0 /)
     ! Increase in Cost of the Audit policy (0.63 is the cost if the pr_audit is 0.5)
     pc_cost_prau_vec = pc_cost_prau_vec - 0.63d0
     ! In the vector below, We start with the cost for the audit prob 82% because the optimal prob is 85% and we need a vector with 18 elements to be consistent with above
     pc_cost_prau_3rd_vec(:) = (/ 0.475640d0, 0.480240d0, 0.487130d0, 0.491730d0, 0.498620d0, 0.503220d0, 0.510110d0, 0.514710d0, 0.521600d0, 0.526190d0, .533090d0, 0.537680d0, 0.544580d0, 0.549170d0, 0.556070d0, 0.560660d0, 0.567560d0, 0.572150d0 /) - 0.02898d0

     ALLOCATE(per_ch_fine_vec(n_chfine))
     per_ch_fine_vec(:) = (/ ( -1d0+0.5d0*DBLE(i) , i=1,n_chfine) /)

     ALLOCATE(prob_conv_vec(n_prconv))
     prob_conv_vec(:) = (/ 0.5d0, 0.76d0, 1d0 /)
     
     ! We set flag_change_ability_policy = 1 if we want to simulate the policies with higher or lower mean ability
     flag_change_ability_policy = 0
     ! ch_mu_a determine the change in mean ability (0.25 or -0.25)
     ch_mu_a = 0d0
     
     ! We change the prob of college degree by +-0.25 to check the robustness of your results; for any other job the change should be zero
     ch_edu = 0d0
     probedu(2) = probedu(2)*(1d0+ch_edu)
     probedu(1) = 1d0 - probedu(2)
     ! We set flag_change_age_policy = 1 if we want to simulate the policies with higher or lower mean agage
     flag_change_age_policy = 0
     ! ch_m_age determine the change in mean age (1 or -1, it corresponds to +- 4 years)
     ch_m_age = -1
     ! We want to increase or decrease the average age by 1.5; in our model 1 period corresponds to 4 years; so we increase or decrease a mayor's age by 1 with probability 1.5/4
     pr_ch_age = 1.5d0/4d0
     
     ntasks = Nboot_tot
     nprocs_per_task = INT(nprocs/ntasks)
     first_proc = 0

     DO iboot = 1, Nboot_tot
        IF (myrank + 1 > nprocs_per_task*(iboot-1) .AND. myrank + 1 <= nprocs_per_task*iboot) nboot = Nboot_min + iboot - 1
     END DO

     IF (myrank == 0) write(*,'(a,5i5)') 'Nboot_tot, Nboot_max, Nboot_min, ntasks, nprocs_per_task', Nboot_tot, Nboot_max, Nboot_min, ntasks, nprocs_per_task
     
     color = nboot
     key = myrank
     CALL MPI_COMM_SPLIT(MPI_COMM_WORLD, color, key, NEW_COMM, ierr)
     CALL MPI_COMM_SIZE(NEW_COMM, nprocs_new, ierr)
     CALL MPI_COMM_RANK(NEW_COMM, my_new_comm_rank, ierr)

     DO iboot = 1, Nboot_tot

        IF (myrank + 1 > nprocs_per_task*(iboot-1) .AND. myrank + 1 <= nprocs_per_task*iboot) THEN
        
           nboot = Nboot_min + iboot - 1
           IF (myrank == nprocs_per_task*(iboot-1)) first_proc = 1

           ALLOCATE(ubp(nparam), lbp(nparam))
           
           OPEN(22, file='inputs.txt')
           READ(22,*) pp
           READ(22,*) ubp
           READ(22,*) lbp
           READ(22,*) factor_param
           CLOSE(22)

           ! We replace the estimated parameters with the ones that corresponds to a particular simulation of the policy
           IF (Nboot_min > 1 .OR. (Nboot_min == 1 .AND. nboot > 1)) pp(:) = param_draws(:,nboot)

           DO i = 1, nparam
              IF (first_proc == 1) write(*,'(a,2i4,2f14.7)') 'nboot, i, pp, factor_param', nboot, i, pp(i), factor_param(i)
           END DO
           
           DO i = 1, nparam
              IF (first_proc == 1) write(*,'(a,2i4,3f20.10)') 'nboot, i, lbp, pp, ubp', nboot, i, lbp(i), pp(i), ubp(i)
           END DO
           
           DO i = 1, nparam
              param(i) = LOG((pp(i)-lbp(i))/(ubp(i)-pp(i)))
           END DO

           ! We need param = pp and flag_sa = 1 to use fcn, if not the parameters are transformed
           param = pp
           flag_sa = 1

           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           !We first run the code for the BASE CASE to get the base past mayor values!
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

           flag_base_policy = 1
           ELNqcons_base_vec_t = 0d0
           ELNqcons_base_vec = 0d0
           policy_pos = Npolicies
           CALL set_policy_vars

           !Read in the matrix of initial conditions
           !the matrix characterizing the nonparameteric distribution
           !the nonparametric distribution for emax
           CALL read_ini_data
           
           ! We run the policies
           CALL fcn(nparam, param, fval)

           flag_base_policy = 0
           DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)
           
           ! We set flag_no_change_ability = 1 if we want to simulate the policies when we draw ability from the unconditional distribution for second term mayors
           flag_no_change_ability = 0
           
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           !We then run the code for the NO-RUN POLICY!
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

           !We then run the code for the policy we are interested in
           flag_norun_policy = 1
           policy_pos = 1
           CALL set_policy_vars
           
           ! We run the policies
           CALL fcn(nparam, param, fval)

           flag_norun_policy = 0
           DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           !We then run the code for the TERMS POLICY!
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           
           !We then run the code for the policy we are interested in
           flag_terms_policy = 1
           policy_pos = 2
           CALL set_policy_vars
           
           ! We run the policies
           CALL fcn(nparam, param, fval)

           flag_terms_policy = 0
           DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           !We then run the code for the TERMS NO-RUN POLICY!
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           
           !We then run the code for the policy we are interested in
           flag_terms_norun_policy = 1
           policy_pos = 3
           CALL set_policy_vars
           
           ! We run the policies
           CALL fcn(nparam, param, fval)

           flag_terms_norun_policy = 0
           DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           !We then run the code for the AUDIT POLICY with cost of audit!
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           
           !We then run the code for the policy we are interested in
           flag_prau_policy = 1
           policy_pos = 4
           cost_policy = 1
           CALL set_policy_vars
           
           ! We run the policies
           CALL fcn(nparam, param, fval)

           flag_prau_policy = 0
           cost_policy = 0
           DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           !We then run the code for the WAGE POLICY with cost of increased wages!
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           
           !We then run the code for the policy we are interested in
           flag_wage_policy = 1
           policy_pos = 5
           cost_policy = 1
           CALL set_policy_vars
           
           ! We run the policies
           CALL fcn(nparam, param, fval)

           flag_wage_policy = 0
           cost_policy = 0
           DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           !THE EVALUATION OF THE REMAINING POLICIES REQUIRE THE COMPUTATION OF THE OPTIMAL AUDIT PROBABILITY FOR DIFFERENT POLICIES DONE EARLIER!
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           ! This is needed to avoid having data written out for the following policies
           flag_optimal_auditprob = 2
           
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           !We then run the code for the OPTIMAL AUDIT POLICY!
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

           ! We compute the optimal audit policy
           flag_optimal_prau_policy = 1
           cost_policy = 1

           ! We computed the optimal audit policy in a previous step
           i = 13
           pr_audit_out = pr_audit_vec(i)
           pc_cost_prau_out = pc_cost_prau_vec(i)
           !We then run the code for the policy we are interested in
           policy_pos = 6
           CALL set_policy_vars
           
           ! We run the policies
           CALL fcn(nparam, param, fval)
           
           DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

           flag_optimal_prau_policy = 0
           cost_policy = 0
           
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           !We then run the code for the TERM LIMITS WITH OPTIMAL AUDIT POLICY!
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

           ! We compute the optimal audit policy with term limits
           flag_term_optimal_prau_policy = 1
           cost_policy = 1

           ! We computed the optimal audit policy in a previous step
           i = 11
           pr_audit_out = pr_audit_vec(i)
           pc_cost_prau_out = pc_cost_prau_vec(i)
           !We then run the code for the policy we are interested in
           policy_pos = 7
           CALL set_policy_vars
           
           ! We run the policies
           CALL fcn(nparam, param, fval)
           
           DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

           flag_term_optimal_prau_policy = 0
           cost_policy = 0

           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           !We then run the code for the NO-RUN WITH OPTIMAL AUDIT POLICY!
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

           ! We compute the optimal audit policy
           flag_norun_optimal_prau_policy = 1
           cost_policy = 1

           i = 7
           pr_audit_out = pr_audit_vec(i)
           pc_cost_prau_out = pc_cost_prau_vec(i)
           !We then run the code for the policy we are interested in
           policy_pos = 8
           CALL set_policy_vars
           
           ! We run the policies
           CALL fcn(nparam, param, fval)
           
           DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

           flag_norun_optimal_prau_policy = 0
           cost_policy = 0

           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           !We then run the code for the TERM AND NO-RUN WITH OPTIMAL AUDIT POLICY!
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

           ! We compute the optimal audit policy
           flag_term_norun_optimal_prau_policy = 1
           cost_policy = 1

           i = 7
           pr_audit_out = pr_audit_vec(i)
           pc_cost_prau_out = pc_cost_prau_vec(i)
           !We then run the code for the policy we are interested in
           policy_pos = 9
           CALL set_policy_vars
           
           ! We run the policies
           CALL fcn(nparam, param, fval)
           
           DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)
           
           flag_term_norun_optimal_prau_policy = 0
           cost_policy = 0

           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           !We then run the code for the TERM AN NO-RUN WITH OPTIMAL AUDIT POLICY ONLY IN THE LAST TERM!
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

           ! We compute the optimal audit policy
           flag_term_norun_optimal_prau_lastterm_policy = 1
           cost_policy = 1

           i = 6
           ! We have 0.8+pr_audit_vec(i) to have as the starting point an audit prob = 83% (the optimal is 85%)
           pr_audit_out = 0.79+pr_audit_vec(i)
           pc_cost_prau_out = pc_cost_prau_3rd_vec(i)
           !We then run the code for the policy we are interested in
           policy_pos = 10
           CALL set_policy_vars
        
           ! We run the policies
           CALL fcn(nparam, param, fval)
           
           DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)
           
           IF (myrank == 0) write(*,'(a)') 'END TERM AND NO-RUN WITH OPTIMAL AUDIT POLICY ONLY IN THE LAST TERM!'
           
           flag_term_norun_optimal_prau_lastterm_policy = 0
           cost_policy = 0           
           
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           !We then run the code for the WAGE POLICY WITH OPTIMAL AUDIT PROBABILITY!
           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

           ! We compute the optimal audit policy
           flag_wage_optimal_prau_policy = 1
           cost_policy = 1

           i = 8
           pr_audit_out = pr_audit_vec(i)
           pc_cost_prau_out = pc_cost_prau_vec(i)
           
           !We then run the code for the policy we are interested in
           policy_pos = 11
           CALL set_policy_vars
           
           ! We run the policies
           CALL fcn(nparam, param, fval)
           
           DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

           flag_wage_optimal_prau_policy = 0
           cost_policy = 0         






           
!!$           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!$           !We then run the code when we change the MEAN OF THE FINE!
!!$           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!$
!!$           ! We compute the optimal audit policy
!!$           flag_hfine_policy = 1
!!$           DO i = 1, n_chfine
!!$
!!$              per_ch_fine_out = per_ch_fine_vec(i)
!!$              !We then run the code for the policy we are interested in
!!$              policy_pos = i
!!$              CALL set_policy_vars
!!$           
!!$              ! We run the policies
!!$              CALL fcn(nparam, param, fval)
!!$
!!$              DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)
!!$
!!$           END DO
!!$
!!$           per_ch_fine = 0d0
!!$           flag_hfine_policy = 0
!!$           
!!$           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!$           !We then run the code for the NO-RUN POLICY WITH A PROBABILITY OF CONVICTION < 1!
!!$           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!$           
!!$           !We then run the code for the policy we are interested in
!!$           flag_probconv_policy = 1
!!$           DO i = 1, n_prconv
!!$
!!$              prob_conv_out = prob_conv_vec(i)
!!$              policy_pos = n_chfine + i
!!$              CALL set_policy_vars
!!$              
!!$              ! We run the policies
!!$              CALL fcn(nparam, param, fval)
!!$              
!!$              DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)
!!$
!!$           END DO
!!$
!!$           flag_probconv_policy = 0

!!$           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!$           !We run the code for the ESTIMATION CASE we only need this to measaure the effect of changes in ability and electoral parameter!
!!$           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!$           
!!$           flag_estimation = 1
!!$           ELNqcons_base_vec_t = 0d0
!!$           ELNqcons_base_vec = 0d0
!!$           ! We set flag_change_ability_01stdev = 1 if we want to determine the effect of increasing mean ability by 0.1 percent of the standard deviation of ability
!!$           ! flag_change_ability_01stdev = 1 changes only how we compute mu_a in moments.f90, we simulate the base case to determine the effect
!!$           flag_change_ability_01stdev = 0
!!$           ! We set flag_change_ability_2ndterm_mayors = 1 if we want to determine the effect of increasing mean ability to the level experienced by 2nd term mayors
!!$           ! flag_change_ability_2ndterm_mayors = 1 changes only how we compute mu_a in moments.f90, we simulate the base case to determine the effect
!!$           flag_change_ability_2ndterm_mayors = 0
!!$           ! If we set to zero: gamma_1=param(8) (general learning); gamma_2=param(11) (diff public consumption); gamma_3=param(9) (clean audit); gamma_4=param(10) (caught stealing)
!!$           !param(10) = 0d0
!!$
!!$           policy_pos = Npolicies
!!$           CALL set_policy_vars
!!$
!!$           !Read in the matrix of initial conditions
!!$           !the matrix characterizing the nonparameteric distribution
!!$           !the nonparametric distribution for emax
!!$           CALL read_ini_data
!!$           
!!$           ! We run the policies
!!$           CALL fcn(nparam, param, fval)
!!$
!!$           flag_base_policy = 0
!!$           DEALLOCATE(betanewcr, betanornew, betasimrun, betapmayor, betapmayortemp1, betapmayornor, beta_voter_i, beta_voter_c)

           
        END IF

     END DO

     IF (Nboot_max == 1) THEN

        bootstrap_vec = bootstrap_vect

        IF (myrank == 0) THEN
           mu_fine = 0.0942693d0 + pp(2)
           DO i = 1, n_chfine
              IF (per_ch_fine_vec(i) > 0d0) THEN
                 mu_finet = mu_fine+ABS(mu_fine*per_ch_fine_vec(i))
              ELSE
                 mu_finet = mu_fine-ABS(mu_fine*per_ch_fine_vec(i))
              END IF
              write(*,'(a,1i4,10f20.10)') 'policy', i, per_ch_fine_vec(i), mu_finet, bootstrap_vec(i,1,:)
           END DO
        END IF
        
!!$        IF (myrank == 0) THEN
!!$           DO i = 1, n_prconv
!!$              write(*,'(a,1i4,9f20.10)') 'policy', i, prob_conv_vec(i), bootstrap_vec(n_chfine+i,1,:)
!!$           END DO
!!$        END IF

!!$     IF (myrank == 0) THEN
!!$        DO i = 1, Npolicies
!!$           write(*,'(a,1i4,10f20.10)') 'policy', i, ch_mu_a, ch_edu, bootstrap_vec(i,1,:)
!!$        END DO
!!$     END IF

!!$        IF (myrank == 0) THEN
!!$           DO i = 1, Npolicies
!!$              write(*,'(a,1i4,8f20.10)') 'policy', i, bootstrap_vec(i,1,:)
!!$           END DO
!!$        END IF
        
     ELSE

        CALL MPI_ALLREDUCE(bootstrap_vect(1,1,1), bootstrap_vec(1,1,1), Npolicies*Nboot_tot*Nvar_sim, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)

        IF (myrank == 0) THEN
           
           OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simparam_base.dat", STATUS='OLD', ACCESS = 'APPEND')
           policy_pos = Npolicies
           DO iboot = 1, Nboot_tot
              IF (Nnterms == 1) WRITE(28,'(6f20.14)') bootstrap_vec(policy_pos,iboot,:)
              IF (Nnterms == 2) WRITE(28,'(7f20.14)') bootstrap_vec(policy_pos,iboot,:)
              IF (Nnterms == 3) WRITE(28,'(8f20.14)') bootstrap_vec(policy_pos,iboot,:)
              write(*,'(a,2i4,8f20.10)') 'base policy', policy_pos, iboot, bootstrap_vec(policy_pos,iboot,:)
           END DO
           CLOSE(28)
           
           IF (Npolicies > 0) THEN
              OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simparam_norun.dat", STATUS='OLD', ACCESS = 'APPEND')
              policy_pos = 1
              DO iboot = 1, Nboot_tot
                 IF (Nnterms == 1) WRITE(28,'(6f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 2) WRITE(28,'(7f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 3) WRITE(28,'(8f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 write(*,'(a,2i4,8f20.10)') 'policy', policy_pos, iboot, bootstrap_vec(policy_pos,iboot,:)
              END DO
              CLOSE(28)
           END IF
           
           IF (Npolicies > 1) THEN
              OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simparam_terms.dat", STATUS='OLD', ACCESS = 'APPEND')
              policy_pos = 2
              DO iboot = 1, Nboot_tot
                 IF (Nnterms == 1) WRITE(28,'(6f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 2) WRITE(28,'(7f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 3) WRITE(28,'(8f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 write(*,'(a,2i4,8f20.10)') 'policy', policy_pos, iboot, bootstrap_vec(policy_pos,iboot,:)
              END DO
              CLOSE(28)
           END IF
           
           IF (Npolicies > 2) THEN
              OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simparam_terms_norun.dat", STATUS='OLD', ACCESS = 'APPEND')
              policy_pos = 3
              DO iboot = 1, Nboot_tot
                 IF (Nnterms == 1) WRITE(28,'(6f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 2) WRITE(28,'(7f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 3) WRITE(28,'(8f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 write(*,'(a,2i4,8f20.10)') 'policy', policy_pos, iboot, bootstrap_vec(policy_pos,iboot,:)
              END DO
              CLOSE(28)
           END IF
           
           IF (Npolicies > 3) THEN
              OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simparam_prau_cost.dat", STATUS='OLD', ACCESS = 'APPEND')
              policy_pos = 4
              DO iboot = 1, Nboot_tot
                 IF (Nnterms == 1) WRITE(28,'(6f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 2) WRITE(28,'(7f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 3) WRITE(28,'(8f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 write(*,'(a,2i4,8f20.10)') 'policy', policy_pos, iboot, bootstrap_vec(policy_pos,iboot,:)
              END DO
              CLOSE(28)
           END IF
           
           IF (Npolicies > 4) THEN
              OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simparam_wages_cost.dat", STATUS='OLD', ACCESS = 'APPEND')
              policy_pos = 5
              DO iboot = 1, Nboot_tot
                 IF (Nnterms == 1) WRITE(28,'(6f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 2) WRITE(28,'(7f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 3) WRITE(28,'(8f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 write(*,'(a,2i4,8f20.10)') 'policy', policy_pos, iboot, bootstrap_vec(policy_pos,iboot,:)
              END DO
              CLOSE(28)
           END IF
           
           IF (Npolicies > 5) THEN
              OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simparam_optimal_prau.dat", STATUS='OLD', ACCESS = 'APPEND')
              policy_pos = 6
              DO iboot = 1, Nboot_tot
                 IF (Nnterms == 1) WRITE(28,'(6f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 2) WRITE(28,'(7f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 3) WRITE(28,'(8f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 write(*,'(a,2i4,8f20.10)') 'policy', policy_pos, iboot, bootstrap_vec(policy_pos,iboot,:)
              END DO
              CLOSE(28)
           END IF
           
           IF (Npolicies > 6) THEN
              OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simparam_term_optimal_prau.dat", STATUS='OLD', ACCESS = 'APPEND')
              policy_pos = 7
              DO iboot = 1, Nboot_tot
                 IF (Nnterms == 1) WRITE(28,'(6f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 2) WRITE(28,'(7f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 3) WRITE(28,'(8f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 write(*,'(a,2i4,8f20.10)') 'policy', policy_pos, iboot, bootstrap_vec(policy_pos,iboot,:)
              END DO
              CLOSE(28)
           END IF
           
           IF (Npolicies > 7) THEN
              OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simparam_norun_optimal_prau.dat", STATUS='OLD', ACCESS = 'APPEND')
              policy_pos = 8
              DO iboot = 1, Nboot_tot
                 IF (Nnterms == 1) WRITE(28,'(6f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 2) WRITE(28,'(7f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 3) WRITE(28,'(8f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 write(*,'(a,2i4,8f20.10)') 'policy', policy_pos, iboot, bootstrap_vec(policy_pos,iboot,:)
              END DO
              CLOSE(28)
           END IF
           
           IF (Npolicies > 8) THEN
              OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simparam_term_norun_optimal_prau.dat", STATUS='OLD', ACCESS = 'APPEND')
              policy_pos = 9
              DO iboot = 1, Nboot_tot
                 IF (Nnterms == 1) WRITE(28,'(6f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 2) WRITE(28,'(7f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 3) WRITE(28,'(8f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 write(*,'(a,2i4,8f20.10)') 'policy', policy_pos, iboot, bootstrap_vec(policy_pos,iboot,:)
              END DO
              CLOSE(28)
           END IF
           
           IF (Npolicies > 9) THEN
              OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simparam_term_norun_optimal_prau_3rd.dat", STATUS='OLD', ACCESS = 'APPEND')
              policy_pos = 10
              DO iboot = 1, Nboot_tot
                 IF (Nnterms == 1) WRITE(28,'(6f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 2) WRITE(28,'(7f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 3) WRITE(28,'(8f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 write(*,'(a,2i4,8f20.10)') 'policy', policy_pos, iboot, bootstrap_vec(policy_pos,iboot,:)
              END DO
              CLOSE(28)
           END IF
           
           IF (Npolicies > 10) THEN
              OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simparam_wage_optimal_prau.dat", STATUS='OLD', ACCESS = 'APPEND')
              policy_pos = 11
              DO iboot = 1, Nboot_tot
                 IF (Nnterms == 1) WRITE(28,'(6f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 2) WRITE(28,'(7f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 IF (Nnterms == 3) WRITE(28,'(8f20.14)') bootstrap_vec(policy_pos,iboot,:)
                 write(*,'(a,2i4,8f20.10)') 'policy', policy_pos, iboot, bootstrap_vec(policy_pos,iboot,:)
              END DO
              CLOSE(28)
           END IF
           
        END IF

     END IF

     DEALLOCATE(bootstrap_vect, bootstrap_vec)
     
  END IF

  CALL CPU_TIME(t1)

  DEALLOCATE(udraw, rnorm, rnormprob, subutil)

  CALL MPI_BARRIER(MPI_COMM_WORLD, ierr)
  CALL MPI_FINALIZE(ierr)

END PROGRAM Main
